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ABSTRACT 


The  stiffness  of  stepped  bars  is  determined  by  calcu- 
lating deflections  due  to  applied  loads  using  the  finite 
element  method.   Nine  combinations  of  step  height  and  fillet 
radius  are  studied  for  each  of  the  following  cases : 

1.  rectangular  cross  section  under  axial  load; 

2.  rectangular  cross  section  under  bending  load; 

3.  circular  cross  section  under  axial  load; 

4-.   circular  cross  section  under  torsional  load. 
A  correction  parameter,  actually  a  fictitious  axial  dis- 
placement of  the  step,  is  derived.   It  makes  possible  ac- 
curate calculation  of  deflections  of  stepped  bars  using 
elementary  formulas . 
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LIST  OF  SYMBOLS 

A     =  cross-sectional  area  (in.2) 

[B]    =  matrix  relating  the  displacement  vector  to  the  strain 
vector 

d      =  smaller  diameter  for  bar  with  circular  cross  section 
(in.),  or  smaller  height  for  bar  with  rectangular 
cross  section  (in.) 

D     =  larger  diameter  for  bar  with  circular  cross  section 
(in.),  or  larger  height  for  bar  with  rectangular 
cross  section  (in.) 

[D]  =  elasticity  matrix 

e  =  superscript ,  denotes  element  contribution 

E  =  Young's  modulus  of  elasticity  (psi) 

[f]  =  element  nodal  force  vector  (lb.) 

[F]  =  global  nodal  force  vector  (lb.) 

G  =  modulus  of  elasticity  in  shear  (psi) 

h  =  height  of  an  element  (in.) 

i  =  nodal  point  number 

I  =  centroidal  moment  of  inertia  (in1!) 

J  =  polar  moment  of  inertia  (ir.1?) 

[k]  =  element  stiffness  matrix 

I  -    axial  length  (in.) 

M  =  bending  moment  (in. lb.) 

N  =  shape  function 

P  =  axial  tension  (lb.) 

r  =  fillet  radius  (in.),  cylindrical  coordinate 

r,z,0  =  cylindrical  coordinates 

T  =  torque  (in. lb.) 

u  =  axial  displacement  (in.) 
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v  =  radial  displacement  (in.) 

x,y,z  =  rectangular  coordinates 

[a]  =  displacement  vector 

y  =  shear  strain  (in. /in.) 

6  =  elongation  (in.) 

<Su  =  virtual  displacement 

6W  =  virtual  work 

A/D  =  stiffness  correction  parameter 

•e  =  longitudinal  strain  (in. /in.) 

6  =  angle  of  twist  (rad.);  cylindrical  coordinate 

k  =  curvature 

v  =  Poisson's  ratio 

C»n  =  normalized  coordinates 

a  =  normal  stress  (psi) 

t  =  shear  stress  (psi) 

4>  =  bending  angle  (rad.) 
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I .   INTRODUCTION 

There  are  many  mechanical  and  structural  applications 
of  bars  and  beams  which  have  abrupt  changes  in  cross-sec- 
tional area.   In  the  design  of  these  members  their  deflec- 
tions under  expected  loads  must  be  considered.   If  accuracy 
is  not  essential,  the  deflections  can  be  determined  by 
relatively  basic  calculations.   These  calculations  involve 
formulas  based  on  the  simple  stress  distributions  found  in 
constant  section  members.   They  do  not  take  into  considera- 
tion the  reduced  stiffness  due  to  stress  flow  irregularities 
in  the  vicinity  of  the  step.   A  stress  concentration  is  said 
to  exist  at  this  location. 

Through  the  numerical  analysis  technique  called  the 
finite  element  method,  the  deflection  of  stepped  bars  or 
beams  can  be  accurately  determined.   This  method  is  sensi- 
tive to  the  effect  of  the  stress  concentration. 

The  purpose  of  this  study  was  to  analyze  the  overall 
stiffness  of  stepped  bars.   An  additional  goal  was  to  de- 
termine a  correction  parameter  which  could  be  used  in  de- 
sign to  compensate  for  the  reduced  stiffness  due  to  the 
stress  concenxration .   Finire  element  analyses  were  performed 
for  stepped  bars  of  various  dimensions  under  axial,  bending, 
and  torsional  loads.-  The  resulting  deflections  were  deter- 
mined, converted  to  a  common  correction  parameter,  and 
compared . 
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The  recent  emergence  of  the  finite  element  method  along 
with  the  digital  computer  has  made  feasible  the  previously 
prohibitive  calculations  necessary  for  stiffness  analysis. 
The  only  prior  study  of  stiffness  of  stepped  bars  that  this 
author  could  identify  was  accomplished  by  F.  P.  Porter  [1]  , 
His  study  of  torsional  vibration  of  irregular  shafts  in- 
cluded stiffness  calculations  based  on  the  mathematical 
theory  of  elasticity  using  an  approximate  graphical 
construction . 


"umbers  in  brackets  refer  to  the  List  of  References, 


p.  80 
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II.   ELEMENTARY  STRESS  RELATIONSHIPS  FOR  THE  THREE  LOADING 
CASES 

A.  BASIC  ASSUMPTIONS 

For  each  of  the  three  loading  cases  (axial,  bending,  and 
torsional),  the  stepped  bar  was  assumed  to  be  perfectly 
elastic,  homogeneous,  isotropic,  and  to  exhibit  a  linear 
stress-strain  relationship.   In  addition,  it  was  assumed 
that  the  bar  was  initially  straight  and  all  loads  were  ap- 
plied smoothly  to  avoid  the  effect  of  impact  loading. 

B.  THE  AXIAL  LOAD  CASE 

A  bar  subjected  to  equal  and  opposite  tensile  loads  will 
deform  according  to  the  equation 

*  ■  1  (1> 

where  6  is  the  elongation  in  length  I ,  P  is  the  axial  ten- 
sion, A  is  the  cross-sectional  area,  and  E  is  the  modulus 
of  elasticity  of  the  bar  material.   This  equation  is  valid 
provided  that  the  loading  is  applied  uniformly  across  both 
ends  of  the  bar  and  that  the  bar  is  of  constant  cross-section, 

For  a  bar  of  non-uniform  cross  section,  the  elongation 
is  estimated  by  dividing  the  bar  into  segments  of  constant 
cross  section  and  adding  the  segment  elongations.   This  pro- 
cedure does  not  take  into  account  the  reduced  stiffness  due 
to  the  stress  concentration  located  at  the  change  of  cross 
section . 
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Figure  1  is  a  longitudinal  section  of  the  bar  studied. 
The  estimated  elongation  of  the  bar  is  therefore 


•  ■  I '  £  ♦  £  * 


(2) 


A,  and  A_  are  the  respective  cross-sectional  areas. 

The  axial  load  case  was  analyzed  for  both  the  plane 
stress  state  and  the  axisymmetric  stress  state.  In  the 
plane  stress  analysis ,  the  bar  was  considered  to  have  unit 
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Figure  1.   Geometry  of  stepped  bar. 

thickness  ,  therefore  the  cross-sectional  area  equaled  the 
height  of  the  bar.   In  the  axisymmetric  stress  analysis  the 
bar  was  considered  to  be  of  circular  cross  section. 


C.   THE  BENDING  CASE 

For  bars  subjected  to  flexure,  or  pure  bending,  a  mea- 
sure of  the  beam  flexibility  is  the  angle  through  which  the 
ends  of  the  bar  rotate.   The  total  angle  can  be  determined 
from  the  equation 


* 


M£ 
EI 


(3) 
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where  <f>  is  the  sum  of  angles  <J>1  and  <j>„  as  described  in  Fig. 
2.   M  is  the  bending  moment  and  I  is  the  moment  of  inertia 
about  the  horizontal  centroidal  axis  of  the  cross-sectional 


Figure  2.   Bending  angles. 

area.   Equation  (3)  is  derived  from  the  following  equation 
which  relates  the  bending  moment  around  the  neutral  axis  of 
an  elastic  beam  to  the  curvature,  k ,  of  the  elastic  curve. 

K        EI  ' 
This  theory  is  valid  provided: 

1.  Plane  sections  remain  plane  during  bending. 

2.  The  oar  is  of  constant  cross  section. 

3.  Bending  is  such  that  shear  strains  are  negligible. 

As  in  the  axial  load  case  ,  the  stepped  bar  can  be  con- 
sidered in  two  sections,  each  of  constant  cross  section. 
The  estimate  for  the  total  bending  angle  is  therefore 


H  9 

a    M  r   1  +   ?  1 
*  "  E  L  T~       T~   J- 


(4) 
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The  bar  analyzed  was  of  rectangular  cross  section  and 
considered  to  have  unit  thickness.   The  bending  moment  was 
applied  by  a  linearly  varying  normal  stress  across  both 
ends  of  the  bar  exerting  tension  on  one  side  of  the  neutral 
surface  and  compression  on  the  other  side  as  shown  in  Figure 
3. 


Figure  3.   Bending  stress  distribution. 


D.   THE  TORSION  CASE 

The  equation  for  angular  deflection  of  a  uniform,  round, 
cross-section  bar  subjected  to  torsion  is 


9  "  JG 


(5) 


G  is  the  angular  deflection,  T  is  the  torque,  J  is  the  polar 
moment  of  inertia  of  the  cross-sectional  area,  and  G  is  the 
modulus  of  elasticity  in  shear.   It  is  assumed  that  plane 
sections  perpendicular  to  the  axis  of  the  bar  remain  plane 
and  that  a  radial  line  remains  straight  when  the  bar  is 
twisted . 

In  this  study,  one  end  of  the  bar  was  held  fixed  and  a 
torque  was  applied  to  the  other  end  such  that  the  shear 
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stress  varied  linearly  from  zero-  at  the  centerline  to  a 
maximum  at  the  outer  radius  as  shown  in  Fig.  4-. 


Figure  4-.   Torsional  stress  applied  to  circular  cross  section. 

As  before,  the  bar  was  considered  in  two  sections  of  con- 
stant cross  section  and  the  total  angle  of  twist  was  the  sum 
of  the  angles  for  each  section.   Therefore  the  estimate  for 
the  entire  bar  is 


T   r  ^1      h_ 

G    Jl    J2 


]  . 


(6) 
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III.   METHOD  OF  STIFFNESS  COMPARISON 


It  is  well  documented  iA&t  the  deflection  of  a  bar  of 
non-uniform  cross  section  is  slightly  greater  than  would  be 
indicated  by  basic  equations  (2),  (4),  and  (6).   This  de- 
creased stiffness  is  caused  by  the  stress  intensification 
at  the  point  of  discontinuity. 

One  way  of  accounting  for  this  decreased  stiffness  is 
to  consider  the  effect  of  the  stress  intensification  on  the 
bar  as  a  whole  to  be  the  same  as  the  effect  of  moving  the 
location  of  the  discontinuity.   Increasing  the  length  of 
the  smaller  portion  of  the  bar  and  decreasing  the  length  of 
the  larger  portion  by  the  same  amount  would  achieve  the  same 
effect  on  its  stiffness . 

This  apparent  change  in  the  location  of  the  bar  discon- 
tinuity was  calculated  and  used  as. a  measure  to  compare 
stiffness  of  bars  of  different  dimensions.   The  size  of 
this  change  was  defined  as  the  variable  "A"  and  is  illus- 
trated in  Fig.  5. 

— HAh- 


1 

p 

k 


i± — +- 


H 


Figure  5.   Geometric  interpretation  of  stiffness  correction 
parameter. 
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To  calculate  A,  equations  (2),  (4),  and  (6)  were  altered  as 
follows : 

I    -A     £9+A 

6  -  lc  -^-  +  -at-]  (7) 

«,  -  A     I   +A 

4>  =  g  C  4—   +  4—  ]  (8) 

T    £,-A      £9  +  A 

In  each  loading  case,  the  deflections  required  to  cal- 
culate 6  ,  <J) ,  and  8  were  determined  by  digital  computer  finite 
element  stress  analysis.   Knowing  the  load,  modulus  of  elas- 
ticity, and  physical  dimensions  of  the  bar,  the  parameter  A 
was  determined.   The  A  parameter  is  common  to  equations  (7), 
(8),  and  (9)  and  therefore  is  a  means  of  comparing  stiffness 
corrections  for  the  three  types  of  loading.   In  addition,  it 
is  independent  of  the  load  applied  because  the  magnitude  of 
the  deflection  and  the  load  applied  are  directly  proportional, 

Since  the  A  parameter  is  a  means  of  correcting  for  the 
deflect:'?",  of  a  stepped  bar  under  load,  it  was  defined  as 
the  "stiffness  correction  parameter." 


IV.   SOFTWARE 

A.   BACKGROUND 

Before  the  stiffness  correction  parameter  could  be  de- 
termined, the  deflections  of  the  bar  under  each  of  the  three 
loading  conditions  had  to  be  determined.   Although  photo- 
elastic  analysis  and  direct  strain  measurement  are  suitable 
for  determining  the  effect  of  stress  concentration,  the  ideal 
method  for  overall  stiffness  evaluation  is  the  finite  ele- 
ment stress  analysis. 

In  very  brief  and  general  terms,  the  finite  element 
method  is  an  advanced  numerical  method  that  has  been  made 
possible  by  the  high  speed  digital  computer.   When  applied 
to  structures  it  can  be  used  to  perform  two  or  three  di- 
mensional stress  analyses.   The  structure  under  study  is 
divided  into  subdivisions  called  finite  elements.   Adjacent 
elements  are  connected  along  the  entire  inter-element 
boundary.   Specific  points  at  regular  intervals  along  the 
element  boundaries  are  designated  and  dexined  as  nodal 
points  or  nodes.   When  a  load  is  applied  to  the  structure, 
the  displacements  of  each  of  these  nodal  points  are  calcu- 
lated and  used  to  determine  the  stresses. 

For  each  of  the  problems ,  maximum  use  of  existing  com- 
puter programs  was  made.   For  the  plane  stress  analysis  of 
•the  axial  and  bending  load  cases,  a  program  titled  PLISOP 
[2]  was  used.   It  was  developed  at  the  Naval  Postgraduate 
School,  Monterey,  California  and  has  been  in  use  for  several 
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years.   For  the  torsional  loading  case,  a  program  titled 
TORT  2  [3]  was  used.   It  was  xLe.velopecL.at  the  University  of 
Wales,  Swansea,  and  had  been  rmstalied  locally  but  had  not 
been  used  or  completely  debugged.   In  order  to  study  the 
fourth  problem,  the  axially  loaded  bar  of  circular  cross 
section,  the  capability  of  PLISOP  had  to  be  expanded  to 
perform  an  axisymmetric  stress  analysis. 

B.   MESH  GENERATION 

The  network  of  lines  describing  the  division  of  a  struc- 
ture into  finite  elements  is  referred  to  as  a  mesh.   Nor- 
mally a  finer  mesh  more  accurately  describes  the  physical 
structure  under  study  and  results  in  a  more  accurate  stress 
analysis.   Unfortunately,  increasing  the  number  of  subdi- 
visions or  elements  greatly  increases  the  computational  ef- 
fort involved  in  the  problem  solution.   As  a  result,  the 
maximum  capabilities  of  most  general  purpose  digital  com- 
puters are  rapidly  reached  as  meshes  are  refined. 

In  many  stress  analyses,  the  area  of  interest  is  only  a 
small  portion  of  the  structure.   In  these  cases  the  mesh 
can  be  very  fine  in  the  area  of  interest  and  much  coarser 
for  the  rest  of  the  structure.   This  minimizes  the  total 
number  of  elements  yet  achieves  accuracy  where  it  is  desired. 
In  this  study,  however,  the  analysis  of  the  structure  as  a 
whole  was  of  interest.   Therefore  it  was  desirable  to  obtain 
a  mesh  as  uniformly  refined  as  possible  within  computer  pro- 
gram limitations.   The  only  practical  means  of  maximizing 
accuracy  yet  minimizing  the  total  number  of  elements  was  to 
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take  advantage  of  symmetry.   This  was  done  by  dividing  the 
bar  in  half  along  the  longitudinal  axis. 

It  was  found  satisfactory  to  take  the  lengths  I,   ,  £_  of 
the  two  bar  segments  equal  to  the  larger  height  dimension 
D.   This  provided  a  nearly  undisturbed  stress  distribution 
at  both  ends  of  the  bar.   Thus  the  stress  intensification 
resulting  from  the  abrupt  change  of  cross  section  is  negli- 
gible at  the  ends  of  the  bar. 

Based  on  these  considerations,  the  stepped  bar  was 
represented  as  shown  in  Fig.  6.   In  each  loading  case  the 
ratio  of  d/D  and  the  fillet  radius  were  varied.   These  vari- 
ations necessitated  minor  changes  in  the  mesh  representation, 
however,  this  general  discussion  will  refer  to  the  arrange- 
ment of  Fig.  6 . 
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Figure  6.   Finite  element  mesh. 

For  this  basic  geometry,  the  mesh  used  consisted  of 
twenty-four  isoparametric  quadrilateral  elements.   The  con- 
cept of  the  isoparametric  element  requires  lengthy  explana- 
tion and  is  not  considered  essential  to  this  development. 
The  reader  is  referred  to  chapter  eight  of  Rcf.  4,  for  the 

background  of  the  isoparametric  element.   Quadrilateral 
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elements  were  used  as  opposed  to  triangular  elements. 
Both  of  the  computer  programs  used  isoparametric  quadri- 
lateral elements.   Such  elements  can  be  linear,  quadratic, 
or  cubic.   In  this  study,  the  only  curved  portion  of  the 
mesh,  the  fillet,  was  adequately  described  by  quadratic 
elements . 

The  quadratic,  isoparametric,  quadrilateral  element  has 
four  corner  nodal  points,  or  nodes,  and  four  mid-side  nodes. 
The  coordinates  of  each  of  these  nodes  are  part  of  the  in- 
put data  for  both  the  PLISOP  and  the  TORT  2  programs.   The 
task  of  providing  this  coordinate  data  as  well  as  data  to 
describe  how  the  elements  fit  together  is  a  considerable 
one.   Fortunately,  another  computer  program  was  developed 
to  perform  just  this  task.   This  program,  titled  PLIMEG  [5], 
generates  a  mesh  as  specified.   It  requires  only  enough  geo- 
metric data  to  describe  the  overall  area  to  be  discretized. 
The  directions  for  use  of  PLIMEG  are  contained  in  the  program. 

The  twenty-four  element  mesh  was  chosen  as  an  adequate 
representation  of  the  region.   A  later  convergence  study, 
which  used  a  finer  mesh,  confirmed  the  adequacy  of  the 
twenty-four  element  mesh  for  this  primarily  comparative 
study . 

Illustrations  of  the  meshes  used  for  the  various  geo- 
metries are  shown  in  Appendix  A. 

.C.   COMPUTER  PROGRAM  EMPLOYED 

1.   The  PLISOP  program  performs  a  plane  stress  or  a 
plane  strain  analysis  by  the  finite  element  method.   Detailed 
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directions  for  its  use  are  contained  in  the  program.   The 
bulk  of  the  input  data'  esaaa-  fee  i^mpsxasad  Ky  the  mesh  generat- 
ing program,  PLIMEG ,  previously  'mKiuronbd.   The  axial  and 
bending  effects  were  achieved  by  the  loads  applied  and  by 
the  boundary  conditions  imposed. 

The  load  vectors  must  be  derived  externally  and  are  ap- 
plied by  specifically  defining  the  load  and  its  direction 
at  particular  nodal  points  in  the  mesh.   Derivation  of  the 
force  vectors  for  both  axial  and  bending  loading  is  shown 
in  Appendix  B,  sections  1  and  2,  respectively. 

The  boundary  conditions  necessary  for  the  axial  loading 
of  the  bar  half-section  shown  in  Fig.  6  consisted  of  con- 
straining- each  node  along  the  bottom  surface  (bar  neutral 
surface)  to  freedom  of  movement  in  the  longitudinal  direc- 
tion only.   In  addition,  one  of  the  corner  nodes  along  this 
surface  was  pinned  to  prevent  longitudinal  motion  as  well. 
For  the  bending  load,  the  nodes  along  the  bottom  surface 
were  free  to  move  in  the  vertical  direction  only.   The  two 
corner  nodes  along  this  surface  were  pinned  to  allow  only 
rotational  displacement  at  these  two  points. 

2.   The  TORT  2  program  performs  a  stress  analysis  of 
circular  sectioned  solids  with  varying  diameter  using  the 
finite  element  method.   The  directions  for  its  use  are 
contained  in  a  separate  booklet  accompanying  the  program. 
The  same  meshes  were  used  to  represent  the  stepped  bar 
except  the  origin  of  the  coordinate  system  had  to  be  at 
the  lower  right  corner  of  the  bar.   The  origin  for  the  mesh 
used  previously  for  the  PLISOP  program  was  at  the  lower 
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left  corner.   A  relatively  simple  computer  routine  can  be 
used  to  renumber  the  nodal  point  data  provided  by  the  PLIMEG 
program,  but  the  connectivity  data  was  revised  by  hand. 

Loading  is  applied  by  specifically  assigning  values  of 
shear  stress  at  particular  nodal  points. 

3.   To  study  a  bar  of  circular  cross  section  under  axial 
load,  an  axisymmetric  stress  analysis  capability  was  re- 
quired.  The  most  direct  way  of  obtaining  this  capability 
was  to  alter  the  PLISOP  program.   Due  to  the  symmetry  of  a 
body  of  revolution,  the  problem  remains  mathematically  two 
dimensional.   The  principal  difference  between  the  plane 
stress  and  the  axisymmetric  stress  problems  is  the  intro- 
duction of  a  fourth  component  of  stress,  the  "hoop"  stress. 
The  details  of  the  development  of  the  axisymmetric  capability 
are  included  in  Appendix  C. 

Since  the  PLISOP  program  is  essentially  unchanged  by 
this  added  capability,  the  directions  for  its  use  are  the 
same.   Appendix  D  contains  notes  on  the  axisymmetric  capa- 
bility intended  to  amplify  the  directions  contained  in  the 
program.   Appendix  E  contains  a  listing  of  the  PLISOP  pro- 
gram in  its  revised  form. 

As  in  the  plane  stress  analyses  of  the  axial  and  bending 
load  cases,  the  load  vectors  must  also  be  derived  externally 
for  the  axisymmetric  stress  analysis.   This  derivation  is 
included  in  Appendix  B,  section  3. 
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V.   RESULTS 

A.  PRESENTATION  OF  RESULTS 

For  each  of  the  loading  cases ,  nine  variations  of  the 
bar  dimensions  were  analyzed.   The  dimensions  and  results 
were  nondimensionalized  by  dividing  each  by  the  dimension, 
D,  which  was  held  constant  throughout  the  study.   The 
values  of  the  stiffness  correction  parameter,  A/D,  for  each 
of  the  four  problems  are  presented  in  tabular  form,  Tables 
I  -  IV,  and  graphically,  Figs.  7-10. 

From  the  graphical  presentations  it  is  evident  that  the 
largest  correction  parameters  are  required  when  there  is  no 
fillet.   As  the  fillet  radius  is  increased,  A/D  decreases. 
In  addition,  as  the  dimension,  d,  is  varied  with  r  held 
constant,  the  largest  correction  is  required  for  the  ratio 
of  d/D  equal  to  one-half. 

The  negative  values  of  A/D  indicate  the  bar  is  stiffer 
than  equations  (2),  (4),  and  (6)  would  show.   This  is  ex- 
plained by  the  fact  that  these  equations,  as  applied,  do 
not  take  into  consideration  the  increased  cross-sectional 
area  resulting  from  the  fillet.   The  fillet  area  was  not 
included  in  the  equations  for  calculating  the  correction  in 
order  to  keep  the  calculations  simple  and  uniform  for  all 
cases . 

B.  COMPARISON  OF  THE  FOUR  CASES 

In  order  to  compare  the  four  cases,  a  plot  of  A/D  versus 
r/D  was  prepared  for  each  of  the  three  ratios  of  d/D.   Figures 
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Table  I 


Table  II 


A/D  for  Axial  Load 
(rectangular  cross  section) 


A/D  for  Axial  Load 
(circular  cross  section) 


\  d/D 
r/D^V 

0.25 

0.50 

0.75 

0.0 

0.19 

0.22 

0.16 

0.12  5 

n.i? 

n.is 

n.nq 

0.250 

0.04 

0.06 

0.0 

\d/D 
r/D\ 

0.25 

0.50 

0.  75 

0.0 

0.10 

0.20 

0.12 

0.125 

0.10 

0.13 

0.05 

0.250 

-0.05 

0.04 

-0.04 

Table  III 

A/D  for  Bending 
(rectangular  cross  section) 


Table  IV 

A/D  for  Torsion 
(circular  cross  section) 


\d/D  !      I 

r/D  \  0.25  ;  0.50 


0.75 


1             | 

0.0          .  0.07    |   0.12 

0.11 

0.125 

0.0 

0.06 

0.04 

0.250 

-0.09 

-0.04 

-0.05 

vd/D 


r/D    \ 

0.25 

0.  50 

0.75 

0  .   U 

0.03 

0.  06 

0.  06 

0.125 

[-0.04 

0.0       ■ 

-0.01 

0.250 

-0.12 

-0.09 

-0.09 
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Figure  7.   Stiffness  correction  for  axial  load  (rectangular 
cross  sect  i  onj . 
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Figure  8.   Stiffness  correction  for  axial  load  (circular 
cross  section) . 
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Figure  9.   Stiffness  correction  for  bending  (rectangular 
cross  section) . 
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Figure  10.   Stiffness  correction  for  torsion  (circular 
cross  section ) . 
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11,  12,  and  13  show  the  effect  on  stiffness  of  the  three 
different  types  of  loading  and  two  types  of  cross  section. 
The  most  severe  type  of  loading  is  clearly  axial  as  indi- 
cated by  the  larger  values  of  A/D.   I.  M.  Allison  [6]  also 
found  the  axial  loading  to  be  the  most  severe  of  the  three. 
His  conclusions,  however,  were  based  on  stress  concentra- 
tion factors. 

Comparison  of  the  two  types  of  cross  section  subjected 
to  axial  load  indicates  the  circular  cross  section  is 
stiffer  than  the  rectangular  cross  section  for  a  particular 
ratio  of  d/D. 

C.   CONVERGENCE  AND  UNCERTAINTY 

The  normal  methods  for  proving  convergence  of  a  finite 
element  solution  are  by  further  refining  the  mesh  or  by 
using  higher  order  elements.   Both  of  these  methods  were 
investigated.   For  the  basic  bar  case  (D=4,  d=2 ,  r=0)  under 
axial  load,  the  convergence  test  results  are  presented  in 
Fig.  14. 

These  results  show  the  improvement  in  accuracy  gained 
by  using  the  quadratic  rather  than  the  linear  elements.   Al- 
though it  can  be  assumed  that  cubic  elements  would  give 
further  improvement,  other  considerations  prevailed.   Since 
it  was  desired  to  cover  a  wide  range  of  problems ,  the  addi- 
tional computational  effort  required  for  cubic  elements  was 
deemed  an  unnecessary  luxury.   Also,  since  the  geometry  of 
the  bar  was  relatively  simple,  consisting  of  plane  or  cylin- 
dric  surfaces,  cubic  elements  would  not  provide  the  same 
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Figure  11.   Stiffness  correction  for  various  loads  with 
d/D  =  .25. 
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Figure  12.   Stiffness  correction  for  various  loads  with 
d/D  =  .50. 
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Figure  13.   Stiffness  correction  for  various  loads  with 
d/D  =  .75. 
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Figure  14.   A/D  vs.l/(no.  of  nodes). 

percentage  improvement  that  the  quadratic  elements  did. 

In  order  to  check  the  accuracy  achieved  by  mesh  refine- 
ment with  quadratic  elements ,  the  PLIS0P  program  limits  were 
altered  to  accommodate  a  ninety-six  quadratic  element  mesh. 
Figure  1L  indicates  this  mesh  has  come  very  close  to  con- 
vergence .   Based  on  these  results ,  the  uncertainty  in  the 
values  of  A/D  in  Tables  I  -  IV  is  considered  to  be  the  dif- 
ference between  the  ninety-six  element  result  and  the  twenty- 
four  element  result,  rounded  to  the  second  decimal  place. 
The  uncertainty  is  therefore  +.02. 

An  additional  test  problem  was  run  to  verify  that  the 
geometry  studied,  with  dimensions  &,  and  £„  equal  to  the 
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larger  height  dimension,  was  adequate  to  ensure  that  the 
stress  intensification  at  the  discontinuity  did  not  signifi- 
cantly affect  the  stress  distribution  at  the  ends  of  the  bar. 
To  check  this  aspect ,  the  length  dimensions  1-.    and  &„  were 
increased  to  1.5  times  the  larger  height  dimension.   The 
resulting  A/D  parameter  was  within  .001  of  the  corresponding 
case  with  the  shorter  length,  thus  confirming  the  sufficiency 
of  the  ratios  Jl  /D=£2/D  =  1.0. 

The  one  geometry  and  loading  case  which  can  be  compared 
with  Porter's  work  agrees  very  favorably.  For  torsion  of  a 
bar  with  d/D  =  0.75  and  r/D  =  0.0,  the  stiffness  correction 
parameters  agree  to  the  third  decimal  place. 
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VI.   CONCLUSIONS  AND  RECOMMENDATIONS 

A.  CONCLUSIONS 

The  stiffness  correction  parameters  obtained  are  useful 
and  applicable  in  the  design  of  stepped  bars.   The  tables 
and  curves  included  in  Section  V  can  be  used  with  confidence 
to  the  first  decimal  place  of  A/D  within  the  range  of  physi- 
cal bar  dimensions  studied. 

The  results  also  yield  important  comparative  data  for 
stepped  bars  of  various  dimensions  under  the  three  types  of 
loading.   The  following  comparative  observations  were  made: 

1.  The  value  of  A/D  for  stepped  bars  subjected  to  axial 
load  is  larger  than  the  corresponding  value  for  bending  or 
torsion . 

2 .  The  stepped  bar  with  ratio' d/D  equal  to  one-half 
yields  a  larger  value  of  A/D  than  the  other  ratios  studied. 

3.  Fillets  reduce  stress  intensification  at  the  loca- 
tion of  the  step  and  therefore  increase  the  stiffness  of 
the  bar. 

In  addition,  it  was  observed  that  the  localized  stress 
intensification  had  a  negligible  effect  on  the  cross-sectional 
stress  distribution  at  a  distance  equal  to  the  height  of  the 
bar.   This  illustrates  the  St.  Venant  principle  of  rapid 
dissipation  of  localized  stresses. 

B.  RECOMMENDATIONS 

It  would  be  of  interest  to  isolate  one  loading  case  and 
through  use  of  more  refined  meshes  and  cubic  elements  achieve 
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correction  parameters  accurate  to  the  third  decimal  place, 
In  addition,  a  great  deal  mor.p  -carapajsative  information 
could  be  gathered  with  a  three  cftmeiiBix>nal  finite  element 
solution  program.   Both  of  these  would  require  altering 
existing  computer  programs  and  would  greatly  increase  the 
computer  core  storage  space  required  for  each  run  in  this 
study. 
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APPENDIX  A:   FINITE  ELEMENT  MESHES  EMPLOYED 
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APPENDIX  B:   DERIVATION  OF  EQUIVALENT  NODAL  FORCES 

1.   Plane  Stress  Analysis  of  the  Axially  Loaded  Bar 

Boundary  stresses  and  distributed  loads  are  converted 
to  equivalent  nodal  forces  by  the  principle  of  virtual  work. 
For  the  case  of  external  stress  a   applied  along  the  boundary 
of  an  element,  the  nodal  forces  are  developed  as  follows: 

6W  =  [<5u.]T[f.]e  =  /  a«6u.d(area) 

substituting 

<5u   =    [6ui]TCNi]T 
[f^6    =    /    a[Nifd(area)  (10) 

where 


6W  =  increment  of  virtual  work 


f .  =  nodal  forces 
1 


a    =    applied  stress 

6u.  =  incremental  nodal  displacements 

6u  =  the  virtual  displacement  of  a  typical  point 
on  the  element  surface 

N.  =  nodal  shape  functions 

e  =  superscript,  denotes  element  contribution. 

The  quadratic,  isoparametric,  quadrilateral  element 
pictured  below  was  used  throughout  this  study.   The  shape 
functions  are  expressed  in  terms  of  the  local,  normalized 
coordinates  (£,n).   For  loading  as  shown,  the  nodal  force 
equations  become : 
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Figure  24.   Quadratic,  isoparametric,  quadrilateral  element. 
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where 

h  =  height  of  the  element 
d(area)  =  dy  =  (h/2)dri  for  bar  of  unit  thickness. 
Substituting  the  shape  functions  and  integrating,  the  result 
for  a    equal  constant  is : 


on 
2 


(11) 


r2 

f 

L  _1' 

Where  nodes  are  shared  by  two  adjacent  elements,  the 
nodal  forces  are  summed. 
2..   Plane  Stress  Analysis  of  a  Bar  Under  Bending  Load 

The  equivalent  nodal  forces  for  the  linearly  varying 
stress  were  derived  in  the  same  manner  as  the  axially  loaded 
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case  except  for  the  applied  stress.   The  linearly  varying 
stress  shown  on  the  element  surface,  shown  below,  can  be 

1=  1 


TJ=  0 


7,=-l 


Figure  25.   Elemental  linearly  varying  stress  distribution, 


represented  by  the  equation 


a  .  ai(  izn  ,  +  ^  1+H  , 


(12) 


Substituting  into  the  equation  for  nodal  forces, 


r  k<  ^ 

-i 


?>    +    <V 


¥f 


■-  1J 


N. 


N, 


N 


dn 


L   1J 


.f.J 


h 
6 


2al    +    2a3 


a. 


(13) 


Applying  this  formula  to  one  end  of  a  particular  bar 
spanned  by  two  elements  as  shown  in  Fig.  26  and  letting  ol-a, 


then  aj 


.2  _ 


a/2  and  o: 


0 


3  "  ""  " 1 

Therefore,  the  nodal  forces  along  this  surface,  num- 
bered from  top  to  bottom  are : 
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3.   Axisymmetric  Stress  Analysis  of  the  Axially  Loaded  Bar 
In  this  case ,  the  region  associated  with  an  element  is 
the  sector  of  an  annulus  with  a  one  radian  central  angle. 
Figure  2  7  illustrates  the  coordinates  used  with  the  axisym- 
metric element  and  the  node  ^,ir":bering  used  in  development 

below . 
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Figure  27.   Axisymmetric  element 
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Beginning  with  equation  (10),  the  equivalent  nodal  forces 
for  a  uniform  external  stress  are  developed  as  follows: 


Cf.]   =  /a»[Ni3d(area) 


(10) 


d(area)  =  r  (1  radian)  dr  =  r  dr 


r  =  r1N1  +  r2N2  +  r3N3 


dr  =  (■r1U[   +  r2N^  +  r^N^dn 


where  N!  =  3N-/3n.   Then, 
l      i 

r3  "  rl 
dr  =  s drl  • 


Therefore , 

d(area)  = 


r3  2  ri  (rlNl  +  r2N2  +  r3N3Mn  * 


Substituting  back  into  equation  (10), 


r3"rl 


/(N1)2dn  /N1N2dn  /N1N3dn' 
/(N2)2dn  /N2N3dn 

Lsym.  /(N3)2dn 


~rl~ 

r2 

Lr3J 

integrating  the  shape  functions  with  g=l,  leaves 


'1 

f2 
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"  4       2 
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r1     2 
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(15) 
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APPENDIX  C:   DEVELOPMENT  OF  THE  AXISYMMETRIC  STRESS  CAPABILITY 

For  an  axisymmetrically  loaded  bar,  displacements  are 
limited  to  the  axial  and  radial  directions.   In  this  study 
these  displacements  were  identified  with  the  letters  u  and 
v,  respectively.   Cylindrical  coordinates  are  applicable 
and  are  defined  in  Fig.  2  8.   It  can  be  seen  that  radial  dis- 
placements introduce  a  circumferential  component  of  strain 


Figure  2  8.   Cylindrical  coordinates. 

as  well' as  a  radial  component.   For  problems  of  this  type, 
the  general  strain-displacement  relationships ,  as  developed 
in  chapter  five  of  Ref.  4,  are 

3v/3r 

v/r 

3v/3z  +  3u/3r 


r  e  ~\ 

z 


[e]  = 


Y 


rz  J 


(16) 


The  elasticity  matrix  which  links  the  strains  to  the  stresses 
is  derived  from  the  following  stress-strain  relationships: 

£z  '-    k    loz    ~   v(or  +  °0)] 

e   =  -  [a   -  v(a   +  a„ ) ] 

r    z,         r       z     0 
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i  [aQ  -  v(ar  ♦  az>] 


Y    =  t   /G 
'rz    rz 


The  resulting  elasticity  matrix  [D]  is 


CD] 


E(l-v) 
(l+v)(l-2v) 


1   v/(l-v)  v/(l-v)    o 

1     v/(l-v)    o 

sym.        1       0 


(l-2v)/2(l-v) 


(17) 


The  finite  element  formulation  of  the  element  stiffness 
matrix  is  carried  out  as  in  the  plane  stress  problem,  ex- 
cept the  integration  over  the  volume  of  the  element  must 
take  into  consideration  the  fact  that  the  element  thickness 
varies  with  radius.   The  element  stiffness  matrix  [k]  can 
therefore  be  expressed  by  the  equation, 


[k]e  =  /  [B]T[D][B]r  dr  dz 


(18) 


where  [B]  is  the  matrix  derived  from  the  strain-displacement 
relationships  and  defined  such  that 

[e]  -    [B][a3 

where  [a]   is  the  element  displacement  vector. 

The  primary  changes  to  the  PLISOP  program  were  to  in- 
corporate the  new  [B],  [D] ,  and  [k]  matrices,  as  defined 
above.   The  element  shape  functions  also  had  to  be  added, 
since  only  their  derivatives  were  included  PLISOP. 

For  further  detail  on  axisymmetric  stress  analysis  by 
the  finite  element  method,  the  reader  is  referred  to  chapter 
five  of  Ref.  4. 
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APPENDIX  D:   NOTES  ON  AXISYMMETRIC  STRESS  ANALYSIS  CAPABILITY 
OF  PLISOP 


1.  z  and  r  denote  respectively  the  axial  and  radial 
coordinates  and  u  and  v  the  corresponding  displacements. 
(In  PLISOP,  z  =  x  and  r  =  y.) 

2.  The  volume  of  material  associated  with  an  element 
has  a  thickness  equal  to  the  arc  subtended  by  a  one  radian 
angle  vice  the  body  of  revolution  usually  discussed  in  most 
texts.   It  turns  out  the  only  difference  is  a  factor  of  2tt. 

3.  Since  the  axisymmetric  analysis  involves  division 
by  r,  there  cannot  be  any  nodal  points  with  coordinate  r=0. 

-  6 

For  nodal  points  along  the  centerline ,  substituting  r=lxlO 
suffices  to  prevent  the  solution  from  "blowing  up,"  although 
the  accuracy  of  stress  values  near  these  nodal  points  is  ad- 
versely affected. 
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